namespace NBody.Analysis
{
    using System;
    using DataFile;
    
    public abstract class Bin
    {
        public delegate double Function(Particle p);
        public enum Spacing { Normal, Logarithmic, EqualMass };
        
        protected int n = 100;
        protected Table data = null;
        protected Spacing spacing;
            	                
        public int NumBins
        {
            get { return n; }
            set { n = value; }
        }
        
        public Table Data
        {
            get { return data; }
        }
        
        public Spacing BinSpacing
        {
            get { return spacing; }
            set { spacing = value; }
        }              
        
    	public abstract void Analyse(NBodySystem s);
        
        protected void radial_bin_normal(NBodySystem s, Function f)
        {
            Console.Error.WriteLine("Binning in radial shells with normal spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
                                                                                
            double deltaR = (rmax - rmin) / (double)n;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);                                                                                
            data = new Table(n, 2, "Radius, Binned quantity");
            data.Clear();
            
            int[] num = new int[n];
            foreach (Particle p in s)
            {
                int pos = (int)((p.Radius - rmin) / deltaR);
                num[pos]++;
                data[pos, 1] += f(p);
            }
            
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = rmin + 0.5 * deltaR + deltaR * i;
                if (num[i] > 0)
                    data[i, 1] /= (double)num[i];
            }
        }
        
        protected void radial_bin_log(NBodySystem s, Function f)
        {
            Console.Error.WriteLine("Binning in radial shells with logarithmic spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
            
            double deltaR = Math.Log10(rmax / rmin) / (double)n;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);             
            data = new Table(n, 2, "Radius, Binned quantity");
            data.Clear();
            
            int[] num = new int[n];
            foreach (Particle p in s)
            {
                int pos = (int)((Math.Log10(p.Radius) - Math.Log10(rmin)) / deltaR);
                num[pos]++;
                data[pos, 1] += f(p);
            }
             
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = Math.Pow(10.0, Math.Log10(rmin) + 0.5 * deltaR + deltaR * i);
                if (num[i] > 0)
                    data[i, 1] /= (double)num[i];
            }
        }
        
        protected void radial_bin_equal_mass(NBodySystem s, Function f)
        {
            Console.Error.WriteLine("Binning in radial shells with equal mass spacing");
            
            int num_per_bin = s.NumParts / n + 1;
            s.Sort();
            
            double rmax = s[s.NumParts - 1].Radius + 1e-10;
            double rmin = s[0].Radius;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);
            
            data = new Table(n, 2, "Radius, Binned quantity");
            data.Clear();
            
            double f_enc = 0.0;
            double n_enc = 0.0;
            double r_last = rmin;
            int c = 0;
            for (int i = 0; i < s.NumParts; i++)
            {
                f_enc += f(s[i]);
                n_enc += 1.0;
                if ((i != 0 && i % num_per_bin == 0) || i == s.NumParts)
                {
                    double r = s[i].Radius;
                    data[c, 0] = r_last + 0.5 * (r - r_last);
                    data[c, 1] = f_enc / n_enc;
                    f_enc = 0.0;
                    n_enc = 0.0;
                    r_last = r;
                    c++;
                }
            }
        }
        
        // returns sigma^2 !
        protected void radial_bin_normal_dispersion(NBodySystem s, Function f)
        {
            Console.Error.WriteLine("Binning dispersion in radial shells with normal spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
            double deltaR = (rmax - rmin) / (double)n;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);
            
            data = new Table(n, 2, "Radius, Dispersion");
            data.Clear();
            
            int[] num = new int[n];
            double[] v = new double[n];
            double[] v2 = new double[n];
            for (int i = 0; i < n; i++)
            {
                num[i] = 0;
                v[i] = v2[i] = 0.0;
            }
            
            foreach (Particle p in s)
            {
                int pos = (int)((p.Radius - rmin) / deltaR);
                double vv = f(p);
                v[pos] += vv;
                v2[pos] += vv * vv;
                num[pos]++;
            }
            
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = rmin + 0.5 * deltaR + deltaR * i;
                double nn = (double)num[i];
                if (nn > 0)
                    data[i, 1] = v2[i] / nn - v[i] * v[i] / nn / nn;
            }
        }
        
        protected void radial_bin_log_dispersion(NBodySystem s, Function f)
        {
            Console.Error.WriteLine("Binning dispersion in radial shells with logarithmic spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
            double deltaR = Math.Log10(rmax / rmin) / (double)n;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);
            
            data = new Table(n, 3, "Radius, Dispersion, error");
            data.Clear();
            
            int[] num = new int[n];
            double[] v = new double[n];
            double[] v2 = new double[n];
            for (int i = 0; i < n; i++)
            {
                num[i] = 0;
                v[i] = v2[i] = 0.0;
            }
            
            foreach (Particle p in s)
            {
                int pos = (int)((Math.Log10(p.Radius) - Math.Log10(rmin)) / deltaR);
                double vv = f(p);
                v[pos] += vv;
                v2[pos] += vv * vv;
                num[pos]++;
            }
            
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = Math.Pow(10.0, Math.Log10(rmin) + 0.5 * deltaR + deltaR * i);
                double nn = (double)num[i];
                if (nn > 0)
                    data[i, 1] = v2[i] / nn - v[i] * v[i] / nn / nn;
            } 
        }
        
        protected void radial_bin_equal_mass_dispersion(NBodySystem s, Function f)
        {
            Console.Error.WriteLine("Binning dispersion in radial shells with equal mass spacing");
            
            int num_per_bin = s.NumParts / n + 1;
            s.Sort();
            
            double rmax = s[s.NumParts - 1].Radius + 1e-10;
            double rmin = s[0].Radius;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);
            
            data = new Table(n, 2, "Radius, Dispersion");
            data.Clear();
            
            double f_enc = 0.0;
            double f2_enc = 0.0;
            double n_enc = 0.0;
            double r_last = rmin;
            int c = 0;
            for (int i = 0; i < s.NumParts; i++)
            {
                double ff = f(s[i]);
                f_enc += ff;
                f2_enc += ff * ff;
                n_enc += 1.0;
                if ((i != 0 && i % num_per_bin == 0) || i == s.NumParts)
                {
                    double r = s[i].Radius;
                    data[c, 0] = r_last + 0.5 * (r - r_last);
                    data[c, 1] = f2_enc / n_enc - f_enc * f_enc / n_enc / n_enc;
                    f_enc = 0.0;
                    f2_enc = 0.0;
                    n_enc = 0.0;
                    r_last = r;
                    c++;
                }
            }
        }
    }
}   		
